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Contrail formation in aircraft wakes using 
large-eddy simulations 

By R. Paoli f, J. Helie $, T. J. Poinsot % and S. Ghosal || 

In this work we analyze the issue of the formation of condensation trails (“contrails”) 
in the near-field of an aircraft wake. The basic configuration consists in an exhaust en- 
gine jet interacting with a wing-tip trailing vortex. The procedure adopted relies on a 
mixed Eulerian/Lagrangian two- phase flow approach; a simple micro-physics model for 
ice growth has been used to couple ice and vapor phases. Large eddy simulations have 
carried out at a realistic flight Reynolds number to evaluate the effects of turbulent mix- 
ing and wake vortex dynamics on ice-growth characteristics and vapor thermodynamic 
properties. 


1. Introduction 

Contrails are ice clouds generated by water in the exhaust from aircraft engines, form- 
ing the common visible white lines in the sky. As assessed in the special report of the 
Intergovernmental Panel on Climate Change (IPCC report (1999)), they have an im- 
portant environmental impact because they artificially increase cloudiness and trigger 
the formation of cirrus clouds, thus altering climate both on local and regional/global 
scales (see figure 1). This was confirmed in a recent climate study (Travis et al (2002)) 
performed in the three days following the 11th of September 2001, when all american 
aircraft were grounded and there were no contrails over the USA. During this period, 
abnormal and significant temperature differences between day and night (namely, daily 
temperature range” or DTR) were observed in the USA. This confirms that contrails are 
of major importance for global climate change by artificial clouds formation. 

Contrails consist of ice crystals which form mainly by condensation of exhaust wa- 
ter vapor over suitable nucleation sites, like soot particles and sulfur aerosols, emitted 
by aircraft engines (Schumann (1996), Karcher & al (1996)). Background atmospheric 
vapor, which eventually adds to the exhaust content, is responsible for the persistence 
of contrails (Gierens (1996)). Contrails form when the air surrounding a nucleation site 
becomes supersaturated with respect to ice (Appleman (1953), Schumann (1996)). For 
certain atmospheric conditions, this may occur somewhere in the jet plume, as the re- 
sult of the increased humidity due to mixing between hot and moist exhaust gases with 
cold and less humid ambient air. In a temperature/vapor-pressure plane (see figure 2), 
assuming that vapor and heat diffuse at the same rate ((Sc = Pr) and that the flow 
is adiabatic, pure mixing can be graphically represented by a straight (“mixing ) line 
which is completely defined by the two states A (ambient air) and B (exhaust gas). In 
figure 2, supersaturation corresponds to the thermodynamic states lying between Si and 
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Figure 1 . Contrails in South-Eastern France sky (http://eol.jsc.nasa.gov). 



Figure 2. Thermodynamic conditions for ice formation: p a (T) is the saturation curve with re- 
spect to ice. State A represents (cold) atmospheric conditions; state B, (hot) exhaust conditions 
(for the sake of representation, we placed it much closer to p s than what it really is). States Si 
and 52 define the range of supersaturation, p w > p 3 (T). 


S 2 , where p w > p g (T). An ice crystal forms when such a condition is locally satisfied 
in the jet plume and, at the same time, a nucleation particle is present. A two-phase 
flow solver, using a Lagrangian solver for nucleation particles, then represents the most 
suitable approach to deal with ice formation in the contrails. 

Formation and persistence of contrails have been studied during last years, mostly 
in geophysical and atmospheric science literature, through in situ measurements and 
numerical simulations with different level of complexity (see IPCC report (1999) and 
references therein). The main intent there was to characterize the general features of 
the phenomenon on time scales up to minutes from the emission time. Other authors 
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investigated in detail the micro-physics of ice formation (the initial composition of the 
jet contrail, the mechanism of water uptake on the nucleation site, homogeneous versus 
heterogeneous nucleation, etc. see Pruppacher & Klett (1997) for a complete review). 

In the present work we have focused on the simulation of contrail formation and early- 
stage evolution in the near-field of an aircraft wake, i.e., up to a few seconds from the 
emission time. The basic configuration consists of an exhaust jet interacting with a wing- 
tip trailing vortex. This represents a complex issue because of the different phenomena 
involved, such as jet and wake-vortex instabilities, transition to turbulence, two-phase 
flow and ice micro-physics. Large-eddy simulations are well suited to deal with all these 
inherently unsteady processes, at realistic flight Reynolds numbers, and have been used in 
this work. The object of the study is two-fold. First, we aim to characterize, (i) the spatial 
distribution of supersaturation around particles, in the jet plume; and (ii) the effects of 
wake vortex. Then, we test a simple micro-physics model for ice-growth to account for 
(iii) the early stage of contrail evolution and (iv) its influence on the thermodynamic 
properties of water vapor. Governing equations for the two-phase flow model are detailed 
in section 2, results are presented in section 3. In section 4 an alternative approach is 
proposed to account for a general (initial) distribution of soot particles. 


2. Governing equations and ice micro-physics model 

Large-eddy simulations of ice formation are carried out through an Eulerian/Lagrangian 
two-phase flow approach. For the gaseous (carrier) phase we solve fully-compressible 
Navier-Stokes equations together with a transport equation for a scalar field, repre- 
senting the exhaust water vapor, Y w . These equations are filtered spatially so th at any 
variable <t>(x) = [p, pu, pv, pw, p E, pY w ] is decomposed into a resolved part <p{x) and a 
non-resolved (or subgrid-scale) part 0"(s), with <j>(x) = 4>{x) + 4>"{x). For compressible 
flows, we use Favre-filtered variables, defined as 4>{x) = <j>(x) -I- <j>'(x), with <p = p<t>! P- 
Using this approach, dimensionless Favre-averaged equations are 


&P d(pUj) 
dt 8xj 

djpui) dfjMjUj) dp _ dojj 

dt dxj dii Re dx 3 dxj ’ 

d{pE) d[(pE + p)Hj} _ 1 dfjjUj doifUj _ 1 c dqj_ _ dQj 

dt + dij ~ Re dij dxj RePr p dxj dxj ’ 

d(pY w ) djpYyUj) = 1 a ( ay,\ 

dt dxj ReScdxj \ dij j dij 


( 2 . 1 ) 

( 2 . 2 ) 

(2.3) 

(2.4) 


The sub-grid scale (SGS) stress tensor <r — —{pu jUj p ujuf), the SGS heat flux Qj 
pC p Tu 3 - pCpfuj and the SGS scalar flux £, = -{pY w Uj - pY w Uj) are modeled through 
the subgrid-scale eddy-viscosity concept 



^ Psgs Cp dQ e psgs dYw 

Qj ~ P^~ dxj ’ ~ Set dxj ’ 


(2.5) 


where 0 = f - ^j~ a kk is the modified temperature (Metais & Lesieur (1992)) while 

turbulent Prandtl and Schmidt numbers are assumed to be constant. In the present 
simulations, we assumed Pr t = 0.3 and Sc t = 0.3. Sub-grid scale viscosity p Bge is obtained 
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through the Filtered Structure Function model (Metais & Lesieur (1992), Ducros & al. 
(1996)), initially developed in spectral space and then transposed into physical space. 
This model was found to be well suited for the simulation of transitional flows, due to 
the property of insuring no SGS viscosity when there is no energy at the cutoff wavelength 
(Ducros & al. (1996)). Vapor/ice mass coupling in the continuity and scalar equations, 
Co is discussed in the next section. 

2.1. Particles treatment and ice-growth model 

Soot-ice particles are tracked through a Lagrangian solver with point force approximation 
(see Boivin al. (1998) for details). Due to their small size (from tens of nanometers to 
a few microns, see Karcher & al. (1996)) their relaxation time is negligible compared to 
the characteristic times of the filtered gas variables. This allows them to be treated as 
tracers. In addition, due to their high number density (varying approximately between 
10 and 10 11 m 3 (Karcher & al. (1996)), we can carry only packets of particles, or 
“numerical particles”, each one containing a large number n tra ns of real soot-ice kernels. 
A numerical particle is defined as the center of mass, x p7 of n tr ans physical particles. In 
the tracer limit, its motion is completely described by 

dx 

— =&(£)$(£ = Sp) ( 2 . 6 ) 

where u is the (filtered) gas velocity. Using filtered quantities instead of exact un-filtered 
ones is equivalent to neglecting sub- grid dispersion, compared to the resolved, large- 
scale dispersion. This assumption is justified because the Reynolds number is high. Gas 
sources are estimated at the numerical particle positions; afterwards they are projected 
on the Eulerian grid, inversely proportional to the mesh-point distance. This numerical 
method can be viewed as a spatial filtering (Boivin & al. (1998), Helie & al. (2002)). 
Exchange between phases is allowed by the code (two-way coupling). Nevertheless, due 
to the tracer limit and solid/gas mass ratio, drag momentum exchange is negligible. 
Temperature cannot be modified by ice growth by more than a few Kelvin (Schumann 
(1996)), so that thermal coupling is also neglected. Therefore, we only consider mass 
exchange, i.e. vapor condensation on the soot particle. The term Co can be neglected in 
(2.1) because of the small amount (order of a few percent) of water vapor in the exhaust 
gases. On the other hand, in (2.4) it accounts for vapor/ice phase exchange and is related 
to the rate of growth of a single ice crystal, through the quasi-stationary model (Karcher 
& al (1996)) 


dr p J yjj 

dt pS p ’ 


P dr 

w = - 22 *(£ = SpKron. -£-pp Sp. 


P=1 


(2.7) 

( 2 . 8 ) 


S p = 4-KTp and p p are the area of the particle (supposed spherical) and its density, while 
J w is the mass flux of water vapor on the particle surface 

J w = 4nr p D G{r p )(Y w - Y t ) (2.9) 

where D = p/(pSc) is the molecular diffusivity of water vaporin air, Y s is the ma«g 
fraction of gaseous water at ice saturation. The dimensionless collision factor G(r p ) is 




FIGURE 3. Basic configuration of a jet/ vortex 
interaction in the near-field, of an aircraft 
wake. 


Figure 4. Sketch of the computational 
domains for the jet and the interaction 
phases. 


given by the semi-theoretical correlation 


G{r v ) 


1 + K n 


AKr 

3a 


■r 


( 2 . 10 ) 


It depends on the Knudsen number K n - \/r p (ratio of the vapor mean free path to the 
particle radius) and the empirical constant a = 0.1, and was found to give good results 
for quasi-isothermal flows and cases with low heat transfer (see Qu & Davis (2001) for 
details) . Saturation conditions are estimated by Sonntag (1994) as 

p s =pX s = exp ( 602 y 28 ? + 29.32707 + 1.0613868 10~ 2 f 

- 1.3198825 KT 3 * 5 f 2 -0.49382577 Inf) (2.11) 


where f is the filtered gas temperature (heat inertia of particles is neglected). Mass and 
molar fractions at saturation conditions sire related by Y a = X S /(X S +(1—X S ) W air /W w ), 
with W a i T = 28.85 kg/Kmole, W w = 18.01 kg/Kmole. 

The numerical code, NTMIX3D, is a three-dimensional, two-phase flow, finite-difference 
solver. For the gas phase, space discretization is performed by a sixth-order compact 
scheme (Lele (1992)). Time integration is performed by means of a three-stage Runge- 
Kutta method, for both the two phases. The solver is fully parallel, using domain de- 
composition. Simulations are performed using 16 processors on CINES 03000 machine. 


3. Results 

This section describes the results of the simulation of ice formation in the near-field 

of an aircraft wake. The basic configuration is an exhaust jet interacting with a (single) 
wing-tip trailing vortex (Fig 3). We adopted a two-stage simulation which consists in 
first simulating a temporally-evolving jet (“jet phase”) and then its interaction with 
the vortex (“interaction phase”). This procedure has been used previously (Paoli & al. 
(2002), Garnier k al. (2002)) and its validity has been recently discussed by some of the 
authors (Paoli k al. (2002)). A sketch of the computational domain for the two phases 
is shown in figure 4. For the jet phase it has dimensions L x = — 16 Tj and L z = 6 Tj 




Figure 5. Joint PDF of temperature and water vapor partial pressure around passive 
particles, during the jet phase (saturation curve is also displayed); left, t = 0.16 s; right 

t = 0.56 s. 


(z being the jet axis and rj the exhaust jet radius) and consists of 161 x 161 x 61 grid 
points; for the interaction phase, the dimensions axe L x = L y = 30 r, and L z = 6r, with 
301 x 301 x 61 grid points. The Reynolds number, based on the radius, r, = 1 m, and the 
exhaust jet velocity, Wj = 60 m/s, is Re j = 3.2 x 10 6 , while the exhaust Mach number is 
Mj = 0.2. Axial velocity, temperature and vapor mass fraction are initialized according 
to a tanh law (for simplicity, all bars and tildes are omitted) 


w o(r) = 5 [ (wj + w a ) - (wj - w a ) F(r)j ; 


F(r) = tanh 


.4 9 \rj r ) 


(3.1) 


where r is the radial distance from the center; Tj/9 — 10; while subscripts j and a 
indicate exhaust jet and ambient air, respectively (see Paoli & al. (2002) for details). 
Initial temperature and water vapor mass fraction follow the same law 


To(r) = | [(7j + T.) - (Tj - T.) F(r)] , (3.2) 

Y* o (0 = 5 [(Y Wi +Y W J- (Y w . - Y w . ) F(r)] . (3.3) 

In the present simulations we assume no coflow, w a = 0, and no background water 
content, Y Wa = 0; the exhaust water content, taken from available data (Garnier & al 
(1997)), is Y Wj = 0.03. Ambient temperature is T a = 225 K and exhaust-to-ambient 
temperature ratio is Tj/T a = 2 which gives Tj = 450 if, a reasonable value for the 
exhaust temperature. The jet is loaded with n p = 250000 (numerical) soot particles with 
the same radius r p = 0.02 /xm. They behave as tracers and each one represents a packet of 
ntrana = 10 6 physical particles. Particle number density is then 2.5 x 10 u , in agreement 
with literature data (Karcher & al (1996)). 

A random-noise perturbation 6w is added to the base flow wo in (3.1) (with Sw max = 
0.005 u/ 0 (r)), to trigger jet instability and transition to turbulence. When the maximum 
jet velocity has decreased to half the initial value, the domain is enlarged and a vortex 
inserted (the relative position is x jv = 5r ; and y jv = -r, , see figure 4), according to the 
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Figure 6. Passive particles case (jet phase). Plane cut of vapor content and distribution of 
supersaturated particles; left, t = 0.56 s; right, t — 0.7 s. 

Lamb-Oseen model (a = 1.4, j3 — 1.2544) 

v e (r) = avA^l-exp -&(£) }: f = Pj (3 ' 4) 

where r c = r, is the vortex core radius and v c = 1.5 Wj(tj v ) is the core velocity (tj v being 
the time at the end of the jet phase simulations; see Paoli & al. (2002)). 

3.1. Passive- particle results 

A first set of simulations was performed with the ice-growth model switched off. The ob- 
ject was to obtain a reference mixing case at high Reynolds numberss typical of aircraft 
wake configurations. It was also useful to analyze the spatial distribution of supersatu- 
rated particles and identify the regions where ice formation is most likely to occur. Basic 
diagnostics consists in analyzing the thermodynamic properties of the exhaust gas during 
mix ing with ambient air. Probability Density Function (PDF) is a convenient diagnostic 
because it provides the additional information on the fraction of particles that super- 
saturate with respect to ice. Figure 5 shows the joint PDF of normalized temperature, 
(T - T a )/{Tj - T a ), and the partial pressure of water vapor, (p w - Pw a )/(Pw j ~ Pu>J, 
around soot particles. The PDF follows a straight line which indicates pure mixing be- 
tween hot jet and cold air. This is a consequence of the assumptions of low Mach number 
and Le = Sc/Pr = 1.. The first assumption implies small pressure fluctuations and ki- 
netic energy negligible compared to internal energy in (2.3). The second implies that the 
diffusion terms in Eqs. (2.3) and (2.4) are the same. Therefore, T and p w are given by 
the same transport equations and evolve along a mixing line (obtained by elimination of 
r in the initial conditions, Eqs. (3.2) and (3.3)) 

Pw = T 1 ( Pwj+Pw. _ T ]+ T a \ 5) 

Pvjj ~ Pu l a — Ta 2 VPtL'j — Pwa Tj —T a ) 

All particles are initially (at t = 0.16 s) placed below the saturation curve p s because they 
are still concentrated inside the hot jet region. Due to the mixing with cold air, particles 
cool, moving along the mixing line, until some of them become supersaturated (crossing 
of pi curve at t = 0.56 s). The spatial distribution of supersaturated particles is given in 
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FIGURE 7. Passive particles distribution during the jet/vortex interaction phase. Dry soot par- 
ticles are represented in black, iced supersaturated particles in white. Total vorticity iso-surface 
identifies the vortex core and the secondary vortical structures due to the interaction with the 
jet; a), t = 0.7 s; 6), t = 1 s; c), t = 1.5 s; d), t = 2 s. 


figure 6, together with a plane cut of water vapor content at two times during the jet 
phase. The figure shows that air first saturates around the particles placed at the edges 
of the jet where the temperature has fallen and there is sufficient vapor to condense. 
The dynamics of the jet/vortex interaction phase are dominated by the entrainment of 
the jet inside the vortex field. Nevertheless, other phenomena occur, as shown in figure 
7: when the jet is close enough to the vortex core, its axial velocity strongly interacts 
with the vortex tangential velocity, causing the formation of three-dimensional structures 
of azimuthal vorticity. These structures progrssively decay (£ = 2 s, see figure 7(d)), 
corresponding to complete entrainment of the jet (Paoli & al (2002)). This mechanism 
of entrainment enhances mixing with external air: therefore, exhaust cooling and vapor 
condensation are favored by the presence of the vortex. Figure 7 shows that at t = 2 s 
all particles are supersaturated and contrail can form everywhere in the wake. 







p w (Pa) 
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Figure 8. Trajectories of three sample Figure 9. Time evolution of three sample 

particles in a T — p w plane. particles radius. 


3.2. Freezing-particle results 

This section presents the results of the simulations with the ice-growth model activated. 
The goal is to analyze the early-stage evolution of the contrail and how it influences 
mixing and the thermodynamic properties of the vapor. The simulation procedure is the 
same as in the previous calculations, i.e., we first simulate the jet alone and then its 
interaction with the vortex. The key point is that the ice and vapor phases are cou- 
pled through Eqs. (2.7) and (2.8). Figure 8 displays the trajectories of three sample ice 
particles in the temperature/ vapor-pressure plane (results are reported only during the 
interaction phase when ice/ vapor coupling is significative). The figure shows that con- 
densation causes large deviations from the mixing line because of vapor removal and the 
consequent decrease in p w . In addition, all the particles finally collapse on the satura- 
tion curve, p s (T), which indicates the thermodynamic equilibrium between vapor and ice 
phases. This is confirmed in figure 9 by the evolution of ice-particle radii which attain 
plateau values between 3 and 6 /im. The distribution of ice crystals size is provided in 
figure 10 in terms of radius PDF. The peak around r v = r p at t = 0.7 s (end of the jet 
phase) indicates that only a small amount of ice has formed. As long as ice nucleation 
proceeds, such a peak decreases and finally disappears, and the shape of the PDF finally 
approaches a Gaussian at approximative^ t = 2 s. A remarkable result is the high vari- 
ance of the PDF with r l p / < r p >« 0.5. This indicates high polydispersion whereas the 
temperature and partial pressure around particles have become approximately uniform. 


4. Future directions 

The particle- tracking approach adopted in this work is an attempt to approximate 
the behavior of a cloud of particles by attributing the properties of a “representative” 
particle to all members of the cluster. In reality, however, every fluid element contains 
particles of many different radii due to different amounts of ice-condensation on each 
of them. Indeed, due to the chaotic nature of particle trajectories in a turbulent flow, 
particles experience quite different ambient conditions during their history. Therefore, 
the appropriate variable that describes the current state of the system is the distribution 
function of particles n p (r, x, t): the number of particles of radius between V and ‘r + dr ’ 
that are contained in am elementary volume ‘dx’ around location ‘x’ at time t. If the law 
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Figure 10. PDF of particles radius during the interaction phase; — o — t t = 0.7s; — o — ? 
t = 1 s; A , t — 1.5 s; - - x - - , t — 2s; ,t = 2s (Gaussian curve). 


of growth of an ice crystal, (2.7) is known, an evolution equation may be written for n v 

&Tlp __ f v 0 , , , 

- + V-(n p u) = -— (n p r). ( 4 . 1 ) 

This is essentially a Liouville equation for particle conservation in the four-dimensional 
phase space ( x,y,z,r ). The term on the right-hand side represents losses or gains from 
each “bin” (r, r+dr ) due to evaporation or condensation. The source term in the equations 
for mass/momentum/energy are of the form / 0 °° (,..)n p dr where (....) denotes the flux of 
the appropriate quantity onto a particle of size ‘r’. In principle, one needs to simulate (4.1) 
for rip, together with the equations of compressible flow. However, since solving partial 
differential equations in four-dimensional space involves a large increase in computational 
cost, it would be prudent to look for simplifying approximations. This is discussed next. 


4.1. Method of moments 

Let us assume that the size distribution of particles at a given location may be written 
as 7i p = F(r;mo,mx,- • •), where mo, mi, ••• sure parameters, which without loss of 
generality, may be taken as the moments of the distribution. The functional form of ‘F’ 
is presumed known. It then remains to determine evolution equations for the moments 
of the distribution, defined as m k {x,t) = r k n p (r,x,t) dr. The first few moments are 
familiar, for example, mo = N p , the total concentration of particles mi = N p (r), where 
(r) is the mean size of particles; m 2 /m 0 - m\/ml = (A r 2 ), the variance of particle size 
distribution about the mean. If we integrate both sides of (4.1) with respect to r, we get 
the following equation for mo = N p 

^ + V.(uN p ) = 0. ( 4 . 2 ) 

This may also be written as D/Dt(N p /p) = 0 using the continuity equation (D/Dt is the 
material derivative). It simply expresses the conservation law for the number of particles 
(irrespective of size). In general, on taking the r-th moment of both sides of (4.1) and 
using integration by parts, we have 


dm k 


rOC 

+ V • (m*u) = / n p r k ~ l f dr. 

Jo 


dt 


(4.3) 
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The right-hand side may be evaluated if a differential equation for particle growth r, and 
the form of the “presumed pdf” n p = F(r ; ■ • ■ ,771*, * * •) is known. If F is assumed to be 
a distribution with ‘n’ parameters, then ‘n’ moment equations need to be retained. 


Example 1 Let us assume the following distribution function; n v = N p 5{r - r p ), that is, 
all particles in a given fluid element at a particular time instant are identical in size. 
Then mi = N p {r) = N p r p , and on substituting this in the moment equation (4.3) we 
have d(N p r p )/dt + V • (7V p r p u) = N p r | r=rp which, using (4.2), may be put in the form 


Dr p 

~Dt 


(4.4) 


Solving (4.4) together with D / Dt(N p / p) = 0 is exactly equivalent to tracking a certain 
number of representative Lagrangian particles released into the fluid (of course, in order 
to obtain the exact solution an infinite number of such particles must be tracked). Particle 
tracking is therefore a special case of our approach. 


Example 2 A more realistic model may be n p = N p <p(r ; mi, m 2 ), that is, the particle size 
distribution at any location is parametrized by its mean and variance. An appropriate 
form for <j> may be chosen on examining the available data in the atmospheric sciences 
literature. Using the growth model (2.7) in (4.3) we then have the following coupled 
equations for determining mi and m 2 

+ V • (m lU ) = —(Y w - r s )Ji(m 1 ,m 2 ), (4.5) 

at pi 

+ V - (m 2 u) = — (Y w - Y.) (4-6) 

at Pi 


where pi is the density of ice and J\, J 2 are functions of mi and m 2 defined by 


7 1 (mi,m 2 ) = / 

J To 


^•d»(r;mi,m 2 )dr, J 2 (mi,m 2 ) = j G(r)4>(r-,m 1 ,m 2 ) dr. (4.7) 
r Jr 0 


Once the form of the presumed pdf and the initial radius of soot particles, ro, is 
selected, the integrals can be evaluated. Clearly, Eqs. (4.2), (4.5) and (4.6) may also be 
put in the Lagrangian form using the continuity equation 


D_ 

Dt 

D_ 

Dt 

D_ 

Dt 



= 0 , 

= - Yv -'~ Ys Ji(m u m 2 ), 
Pi P 

_ D Y W zlL j 2 (m u m 2 ). 
Pi P 


(4.8) 

(4-9) 

(4-10) 


In these equations, u is the velocity vector and all variables denote the physical field. In 
order to obtain the filtered fields, one should apply Favre averages on both sides of Eqs. 
(4.8) - (4.10) and introduce appropriate subgrid modeling assumptions (such as neglect 
of the fluctuating terms and introduction of “eddy” diffusivities) to achieve closure. 


5. Conclusions 

In this work, we studied the process of formation and early evolution of a contrail in 
the near-field of an aircraft wake. The basic tool was a two-phase flow code; the basic 
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configuration, an engine exhaust jet, loaded with soot particles, and interacting with a 
wing-tip trailing vortex. Large-eddy simulations were carried out at high Reynolds num- 
bers, typical of wake configurations. A simple micro-physics model was used to account 
for nucleation of water vapor on the soot particles. A first set of simulations was carried 
out without any ice-growth model, to analyze the spatial distribution of supersaturation 
around soot particles. Results showed that particles first saturate at the edges of the 
exhaust jet. Vortex-induced entrainment in the wake enhances mixing of exhaust gases 
with cold air and favor water vapor supersaturation and condensation. A second set of 
simulations was carried out with the ice-growth model activated. The results showed that 
the radii of the frozen particles reach asymptotic values which depend on the local super- 
saturation. This corresponds to local thermodynamic equilibrium state between vapor 
and ice phases. Taking vapor depletion into account results in significant deviation from 
the classical mixing line. This justifies the need for two-phase flow simulations to deal 
with contrail formation. 
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